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BODIES OF REVOLUTION HAVING MINIMUM 


TOTAL DRAG IN HYPERSONIC FLOW 

By Waldo I. Oehman and Sylvia A. Wallace 
Langley Research Center 

SUMMARY 

The calculus of variations was used to study minimum drag bodies of revolution in 
hypersonic flow. Newton’s formula for pressure drag was assumed, and viscous drag 
was formulated by use of a variable skin-friction coefficient. Constraints of either given 
length and base height or given length and volume were imposed on the bodies. 

Fineness ratios to about 9 were investigated. The minimum drag bodies obtained 
were characterized by flat noses having, at most, a diameter of 3.5 percent of the maxi- 
mum body diameter. The effect of viscous drag was a reduction of the volume of the 
bodies for the larger fineness ratios. Practical body shapes, for which the body slope 
was required to be greater than zero, have an upper limit on the value of fineness ratio. 

By use of a '’slender body'’ approximation, other investigators have obtained a similar 
limit on the fineness ratio. For a few cases investigated for which the ratio of volume 
to body length was held constant rather than the fineness ratio, it was found that viscous 
drag had almost a negligible effect on the minimum drag body shapes. 

INTRODUCTION 

The meridian shapes of bodies of revolution that produce minimum Newtonian pres- 
sure drag have been calculated by Eggers, and others. (See ref. 1.) Suddath and Oehman 
(refs. 2 and 3) calculated the minimum pressure drag shapes for bodies with elliptical 
cross sections. Several authors have recently included the effect of viscous forces in the 
problem of determining the shapes of minimum drag bodies in hypersonic flow. In refer- 
ence 4, Kennet considered slender bodies (that is, bodies for which the local slope is much 
smaller than unity) having an assumed constant skin-friction drag coefficient. Bryson, 
in chapter 18 of reference 6, dropped the slender-body approximation, but retained the 
assumption of a constant skin-friction coefficient. A study by Miele and Cole (refs. 5 
and 6) of two-dimensional shapes and slender pointed bodies of revolution included a 
variable skin-friction coefficient. 

The present investigation extends the previous work by considering nonslender 
bodies, both pointed and flat-nosed, and a variable skin-friction coefficient. In this 


! 


extension a calculus of variations solution of the problem has been obtained, numerically, 
on a digital computer. Computed minimum-total-drag body shapes are presented to 
assess the effect of friction drag and to illustrate the use of the computer program. In 
addition to the constraint of given length and base height, body shapes were computed for 
the constraint of given length and volume. Skin-friction coefficients for either laminar 
or turbulent boundary layers have been used. 

SYMBOLS 


A 


C 


D 


C DF 


c f 


D b 

D f 

D P 

^ '’total 

f,F,g 


constant in skin-friction- coefficient formula 
drag coefficient, — ^ ra f 0 

friction drag coefficient 

local skin-friction coefficient 

base drag of body 

friction drag 

pressure drag 

sum of Df,, Df, and Dp 

integrand functions 


f(x,y,y') integrand function in equation (6) 
F(x,y,y',u,A,j3) integrand function in equation (10) 


F = 

nn 


d*F 


uu " _ o 
du^ 


F » * 

y y 


9 2 F(x,y,y',u,X,ff) 

9 y ' 2 


F < = 

y 


BFCx^y ,u,X,/3) 
9y' 


g(y) integrand function (see eq. (7)) 
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G 


value of integral in equation (7) 


I functional defined by equation (6) 

J functional defined by equation (9) 

L length of body 

n fineness ratio defined as ratio of length to base diameter of body of revolution 

q free-stream dynamic pressure 

OO 

u,y f body slope, ^ 

V volume of body of revolution 

x,y body coordinates, nondimensional with respect to L, along and perpendicular, 

respectively, to axis of rotation 

a exponent in skin-friction-coefficient formula 

j8 Lagrange multiplier 

6 variation consistent with prescribed boundary conditions 

6 variation taken at constant station x 


o constant in skin-friction-coefficient formula 


A 

*(y 0 ) 


multiplier function 

function of body radius at x = 0 (see eq. (6)) 


Subscripts: 


o value at body nose, x = 0 


1 value at body base, x = 1 


i ith value 

Dot over a symbol denotes derivative with respect to x. A prime also denotes a 
derivative with respect to x. 
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PROBLEM FORMULATION AND SOLUTION 


Statement of the Problem 

A body of revolution, at zero angle of attack, is considered to be in a hypersonic 
flow of air. The total drag of the body may be expressed as the sum of pressure drag, 
skin-friction drag, and base drag, or 

D tota! = D p + % + D b <» 


The pressure drag is assumed to satisfy Newton's law of resistance and is formulated as 
follows : 



Formulation of the skin-friction drag is the following integral: 

D f = 27 r( ico L2 J 0 y c f to 

The local skin-friction coefficient is assumed to satisfy 

A 


c f = 



(3) 

(4) 


where the value of the parameter A depends on the physics of the flow conditions, and 
the value of the exponent a depends on the character of the boundary layer. Usually, 
a = 0.2 is associated with a turbulent boundary layer and a = 0.5 is associated with a 
laminar boundary layer. The quantity in the parentheses represents, approximately, the 
distance along the body meridian. The constant o is included to avoid the problem 
associated with an infinite skin-friction coefficient when x = 0. 

Base drag is assumed to be negligibly small and is set equal to zero. Also, the 
effect of boundary- layer thickness on the pressure distributions has not been included in 
equation (2). Furthermore, the boundary layer does not have a transition region. 

Combining equations (1), (2), and (3), the total drag is 


D total = 2 ^~ L2 Jyn 2 + 


'1 r 


2y" 


1 +y’‘ 


(o + x \J\ + y’j 
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or, in nondimensional form, 


D total 2 
2 ^qco L2 Y ° 


+ 



2y 


-3 


1 + y ' 2 (o+xjL+y*f_ 


dx 


(5) 


If the function y(x) is known, the drag factor may be evaluated by per- 

forming the integration indicated by equation (5). However, the problem of interest is to 
determine the function y(x) that makes the drag factor a minimum. In the development 
of equation (5), it has been assumed that the length of the body is known. Therefore, the 
minimization of the drag factor is to be accomplished for a given length and either a given 
base height or a given volume. 


Method of Solution 

The solution of the problem is obtained by applying the calculus of variations and by 
using a digital computer to solve the resulting two-point boundary- value problem. A brief 
development of conditions that the solution must satisfy is presented in functional form. 

For convenience, a quantity I is defined by 

I £ Ptotal 2 = + J f(x,y,y') dx (6) 

2 ^oo L 0 


where 


and 


$ 



2 

o 


f(x,y,y') = y 


2y 


,3 


i+ *' 2 („ + x£77*r 


If the volume of the body is specified, the functional I is to be minimized so that an 
integral of the form 


G = 



dx 


(7) 


has the prescribed value. Furthermore, it will be convenient to substitute u for y f in 
f(x,y,y T ) and introduce the differential equation 

u - y' = 0 (8) 
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as a subsidiary condition. A new functional J which is to be minimized is given by 


J = *(y 0 ) + H(x,y,\a) + A(x) (u - y’) + jSg(y)[ dx (9) 

or 

J = $(y 0 ) + F(x,y,y',u,X,0) dx (10) 

where A(x) is a variable multiplier and /3 is a constant multiplier. Minimization of 
the functional J is equivalent to minimization of I subject to the integral constraint 
(eq. (7)) and the subsidiary condition (eq. (8)). 

Calculating the first variation of J and setting it equal to zero leads to 


SJ = 



X=1 

^1 

r— — 

\9y ay'/ y J 

+ J 

x=0 ^ 

0 

7|£ - |F W + |Z1 g u + 3F g x 

\3y dx &y'J 8u 9x 


dx = 0 


(ID 


where the symbol 6 denotes a variation consistent with prescribed boundary conditions 
and 6 denotes a variation taken at a constant station x. From familiar arguments con- 
cerning the arbitrary variations 6y, 6y, 6u, and 6A, equation (11) leads to the following 

three necessary conditions: 

(1) Satisfaction of Euler equations 

(2) Satisfaction of end conditions 

(3) Satisfaction of Legendre condition. 

Euler equations .- The Euler equations which must be satisfied over the interval 
0 ^ x = 1 are when the length and base height are given: 


dX _ 9F _ 9f 

dx 9y 9y 


dy = 
dx 


u 




8F 

au 




(12a) 
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and when the length and volume are given: 


dx 9y 9y 9 y 



dx 




8F 

au 


-fUx-o 

9u 




(12b) 


It is possible to obtain a first integral of the Euler equations that contains an additional 
variable multiplier. However, since the integration provides no additional information 
for these problems, the integral has been omitted. 

End conditions.- The initial condition that must be satisfied is 


8 $ 


% 


9H =0 

9y' 

J x=0 


(13) 


which is applicable for both problems, and the final condition is 

Fyil = 0 (14) 

-4=i 

which is applicable only for the problem with given length and volume. 

Legendre condition.- The solution y(x) must be obtained so that the Legendre con- 
dition is satisfied in the interval 0 1x^1; that is, 

F U U S 0 < 15 > 

This set of necessary conditions must be satisfied by the functions y(x), u(x), and 
A(x) in order to extremize the functional of equation (9). Furthermore, the resulting 
extremal is called a weak extremal. It has been implied, in this development, that the 
body slope is greater than zero for all values of x. 


Specific Solutions 

Given length and base height .- The integrand function for the problem in which the 
length and base height are given (see eqs. (9) and (10)) is written explicitly as 
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F = y 


(16) 



The Euler equations are 


dA _ _ 2u 3 A 

** 1 + u2 + 




dy 

— = u 
dx 


(17) 


2u 3 (3 + u 3 ) 


aAux 


2 Q/ + 1 

(l + u 3 ) /l + u2 (a + X(/l + u 3 ) 


Applying the end condition at x = 0 (eq. (13)) gives 

*o = - 3 y 0 


!+ A = 0 


(18) 


which is the initial condition for the differential equation for A. The value of y Q is not 
specified and must be chosen so that the integration of equation (18) from x = 0 to x = 1 
yields the given value of y^. Choices of y Q = 0 or y Q > 0 lead to separate conse- 
quences. If y Q > 0, then, by using equations (17) and (18), the value of u Q is unity. If 
y D = 0, the last equation of the system (eqs. (17)) is an identity and does not provide any 
information for evaluation of u Q . However, if the Euler equation of the drag integral 
(eq. (6)) is evaluated at x = 0, an equation in y 0 ' is obtained. The Euler equation is 

_d_ _9f_ _ af_ _ o 
dx 3y» 3y 

where the function f is defined for equation (6). Since f is of the form 

f(x,y,y') = yg(x,y’) 


the Euler equation becomes 




g = 0 


or 


d / 8 g) 
dx \9y' 


+ y 


9 g 

3y' 


g = o 
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When x = 0, y Q = 0 and the first term of the Euler equation is zero. Thus, at x = 0, 


y’ 


9y- 



= 0 


Explicitly, the last form of the Euler equation becomes 

1 2 y 0 ’ 2 ( 3 + y 0 ’ 2 ) 2 y 0 ' 3 ( 1 + y 0 ' 2 ) 


y o 2 

h’o*) 


(‘♦y.'*r 



or more simply, 



which, on expansion, becomes a fourth-degree equation in y Q ' as follows: 


,4 


4a 

A 


.a 


y 0 ’ 3 + 2 y 0 ’ 2 + i = o 


(19) 


Thus, for y Q = 0, this form of the Euler equation evaluated at x = 0 (eq. (19)) provides 
four values for the initial body slope y ’ = u Q . Application of simple algebraic theorems 
shows that two of the roots of equation (19) are real and positive and that the other roots 
are a conjugate complex pair. Only the real roots are significant for the present 
investigation. 

The Legendre condition (inequality (15)) which must hold for all x is 

(20) 


The conditions of the problem require that y(x) > 0 and that u(x) > 0 for 
0 < x £ 1. Consequently, inequality (20) may be solved numerically to obtain u(x) for 
0 < x < 1. When x = 0, y Q must be greater than or equal to zero (y 0 = 0) and u Q 
must be greater than zero (u 0 > 0). Obviously, when y Q = 0, F uu = 0 and any positive 
value of u D satisfies the Legendre condition; and, when y G > 0 any 0 < u Q § \/3 satis- 
fies the Legendre condition. A plot of u as a function of x is presented in figure 1 to 
illustrate the permissible values of u that satisfy the Legendre condition. Thus, any 
body shape, with y(x) > 0, generated from the other necessary conditions will minimize 
the drag integral only if the body slope u is on or between the boundary curves ( F uu - °) 


F 

uu 


y 

jiu[3 - u*) 

1 + u 21 

t(l +u 2 ) 2 


aA 


+ a ) a 2 . i] . 


crx 




S 0 


+ u* 
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of figure 1. If y Q = 0, the slope u Q may take on any positive value, but for 0 < x S 1, 
the slope u must be on or between the boundary curves. 

Given length and volume .- For the problem in which the length and volume are 
given, the volume is the integral of equation (7) with g(y) = y2 or 


G = 



y2 dx 


( 21 ) 


The explicit integrand function of equation (10) is 


F = y 


2u 3 


1 + u^ 



( 22 ) 


The Euler equations are: 

dx _ -2u 3 


dx l + u 2 


(a + xjl + u 3 j 


- 2/3y 


dy 

~ = u 
dx 


2u 2 (3 + u 2 ) 
2 


ckAux 


(i + u2 ) + 


+ A ~ 0 


(23) 


The end condition at x = 0 is the same as equation (18), that is, A 0 = -2y Q ; and the end 
condition at x = 1 (see eq. (14)) requires that 


Xj = 0 


(24) 


Since y 0 and /3 are not specified, they must be chosen so that G has the specified 
value and Xj = 0 when equations (23) are integrated from x = 0 to x = 1. In order 
for X to go to zero at x = 1, /? must be a negative number. The discussion of the 
Legendre condition in the preceding section also applies to this problem. As before, 
u Q = 1 for y Q > 0 and u Q > 0 for y Q = 0. At x = 1, however, the last equation of 
the system (eqs. (23)) may be solved for uj. That is, uj must be a solution of 


yi u i 


^(3 + uj2j 

(l + u x 2 f 



(25) 
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The solution u-^ = 0 is not admissible since the Legendre condition requires uj > 0. 
Consequently, the quantity in parentheses in equation (25) must be zero; and further anal- 
ysis shows that uj must be near zero. Furthermore, the computed value of uj must 
be greater than the value on the lower Legendre boundary at x = 1 in figure 1 for an 
acceptable solution of the problem. 

RESULTS AND COMPUTATIONS 

The preceding development of the necessary conditions from the first variation of 
the drag integral has led to the problem of solving sets of ordinary, first-order differen- 
tial equations (Euler eqs. (17) and (23)). Some, but not all, of the initial and final condi- 
tions are known. Therefore, a two-point bound ary -value problem must be solved. Usu- 
ally, an iterative procedure with some sort of convergence is applied to obtain specific 
solutions. The approach taken in the present investigation has been to consider the 
unknown initial conditions y Q or u Q and 0 as parameters and to obtain solutions of 
the Euler equations for particular values of these parameters. In addition to initial con- 
ditions, the constants A, a, and cr are parameters that reflect the flow conditions and 
a solution (that is, a body shape) may be obtained for each set of parameters. Because 
the numerical results presented in the following sections were obtained for a few selected 
sets of these parameters, the digital computer program is presented, for convenience, in 
the appendix. The computer time is about 3 minutes for each set of initial conditions. 

Given Length and Base Height 

A family of flat-nosed minimum-drag bodies has been computed for the parameter 
sets (A, a, cr) of (0, -, -), (0.0002, 0.5, 0.01), (0.0005, 0.2, 0.01), and (0.0009, 0.2, 0.01), 
and for initial body ordinates y 0 ranging from 10"^ to 10"30. The initial body slope 
was, of course, u Q = 1 for each body of the family. A plot of y 0 as a function of fine- 
ness ratio (that is, the ratio of length to base diameter) is presented in figure 2. The 
minimum pressure drag bodies (A = 0) are characterized by the continued increase of 
fineness ratio as y D tends to zero. However, the fineness ratio of minimum drag bodies 
(A > 0) approaches an upper limit as y 0 tends to zero. The upper limit of the fineness 
ratio for each set (A, a, a) is the value of n indicated for y 0 = 10 - ^ in figure 2. These 
values are also the same as were obtained for y Q = 10" 30, The main result for flat-nose 
bodies that produce minimum drag is that the choice of fineness ratio for given flow con- 
ditions is limited. 

For pointed bodies (y Q = 0), solution of equation (19) yields initial values of the body 
slope y 0 \ The real values, computed for the present investigation, are 
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A, a, a 

V = u o 

0.0002, 0.5, 0.01 
.0005, 0.2, .01 
.0009, 0.2, .01 

0.0797058575; 1999.99900 
.0681782789; 3184.85674 
.0830586184; 1769.36407 


The values of y Q ' depend on the set (A, a, a) which embodies the flight conditions. Fur- 
thermore, since y 0 (y 0 = 0), y 0 '(y 0 ' = u Q ), and A 0 (A 0 = 0) are known, the minimum-drag 
body shape is completely determined by integrating equations (17) from x = 0 to x = 1. 
Consequently, for a given length, the base height cannot be specified arbitrarily but is 
dependent on the flight conditions. For each set (A, a, a), there are two values of y 0 ’ 
each of which leads to a solution (body shape) of the problem. Both body shapes are 
almost identical for x > 10"3. However, the total drag of the body having the smaller 
initial slope is slightly less than the total drag of the body having the larger initial slope. 
Therefore, when solving equation (19) for y Q ', only the smaller value yields the correct 
minimum-drag body slope. It should be emphasized that the restriction on the choice of 
fineness ratio noted here and in the preceding paragraph is a consequence of requiring 
that u(x) > 0. However, the body shapes for larger fineness ratios obtainable when 
u(x) = 0 for some 0 < x < 1 are not compatible with the mathematical model of the 
drag given herein. The restriction on the choice of fineness ratio is the same as may be 
derived from the results given in reference 5 for "slender" bodies with a subcritical value 
of the friction parameters. However, the minimum-drag bodies of reference 5 are blunt 
bodies that have an infinite slope at the nose. 


The body shapes presented in figure 3 for fineness ratios of 2 and 5 illustrate the 
effect of the friction drag. In figure 3(a), the body shapes (n = 2) are almost identical for 
the three sets of (A, a, a) that is, (0, -, -), (0.0005, 0.2, 0.01), and (0.0009, 0.2, 0.01). 

Thus, the body shaping is dominated by the pressure drag. However, for n = 5, 
increasing the local skin-friction coefficient (that is, increasing A) causes a decrease of 
the local body radius. (See fig. 3(b).) Thus, the viscous term in the drag integral leads 
to a decrease of the volume of a minimum drag body with a given fineness ratio. The 
curves in figure 4 show that the friction drag coefficient becomes larger than the pressure 
drag coefficient as the fineness ratio increases. Thus, the importance of including vis- 
cous drag in the problem formulation is emphasized. 


Given Length and Volume 

Sets (A, a, d) of (0.0005, 0.2, 0.01) and (0.0009, 0.2, 0.01) were combined with values 
of y Q ranging from 10“ 2 to 10 - ? (flat-nosed bodies) to compute minimum-drag bodies 
having given length and volume. The calculations were performed by iterating to obtain 
a value of 0 for each set of (A, a, c r) and y Q that would force A to go to zero when 
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x = 1. Figure 5 presents y 0 and the corresponding multiplier /3 for the volume ratio 
V/L 3 . Two points for A = 0 are also presented in figure 5(a). (These points for flat- 
nose bodies are taken from ref. 1.) This figure shows that, for a specific volume ratio, 
a slightly different value of y Q is required for each set (A, a, o). Figure 6 shows that 
the resulting fineness ratio is not appreciably affected by changes of A. Consequently, 
the effect of viscous drag on the body shape is almost negligible although the viscous 
drag coefficient is as much as 80 percent of the total drag coefficient. (See fig. 7.) 

One effect of viscous drag was the requirement imposed by the Legendre condition 
that the body slope at the base must be greater than zero. (Minimum pressure drag bodies 
of refs. 2 and 3 with given length and volume have a zero slope at the base.) However, the 
required slope is very small compared with unity. Figure 8 presents some minimum-drag 
body shapes for given length and volume for several values of the volume ratio. 

CONCLUDING REMARKS 

An analytical investigation was made to determine the meridian shapes of minimum- 
drag bodies having either given length and base height or given length and volume. The 
flow was assumed to be hypersonic and Newton's formula for pressure drag to be appli- 
cable. Viscous drag was formulated by using a skin-friction coefficient that varied with 
distance along the body meridian and the boundary- layer flow was considered to be either 
laminar or turbulent. A solution was obtained by application of the calculus of variations, 
and equations obtained from the necessary conditions were programed for a digital 
computer. 

A limited parametric analysis has been made to assess the effect of viscous drag on 
minimum-drag body shapes and to illustrate the use of the computer program. Minimum- 
drag body shapes, with fineness ratios to about 9, were characterized by flat noses having, 
at most, a diameter of 3.5 percent of the maximum body diameter. The effect of viscous 
drag was a reduction of the volume of these bodies for the larger fineness ratios. Practi- 
cal body shapes, for which the body slope was required to be greater than zero, have an 
upper limit on the value of fineness ratio. A similar limit on the fineness ratio has been 
obtained by other investigators who have used a "slender body" approximation. The main 
conclusion for flat-nose minimum-drag bodies having given length and volume is that the 
viscous drag has almost a negligible effect on the body shapes. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., June 17, 1969, 

126-13-02-27-23. 
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APPENDIX 


DESCRIPTION OF PROGRAM 

The computer program was written in the FORTRAN IV language under SCOPE 
Version 3.0 for Control Data Corporation's Series 6000 computers at the Langley Research 
Center. A brief description of the program, as well as a flow diagram and actual listing 
of the program is included in this appendix. The output from an example problem is also 
given. 


Contents of the Program 

The program numerically integrates the Euler equations (eqs. (23)) which must be 
satisfied over the interval 0 ^ x g 1. Subroutines FALG, INT1A, and ITR2 are used for 
the solution of the quartic (eq. (19)) for the two real values of u Q , for the integration of 
the differential equations, and for the implicit solution of the u (eqs. (17) or (23)), 
respectively. The program determines the drag given in equation (5) using the trape- 
zoidal rule to compute the separate integrals for the pressure drag and friction drag 
coefficients. 


Subprograms 

The program makes use of three library subroutines, FALG, INT1A, and ITR2, from 
the Langley Research Center. A complete listing of these subroutines is included in this 
paper. 

FALG is a subroutine which calculates the n roots of a polynomial of degree n, 
where the coefficients may be either real or complex. In this program it is used to obtain 
the u 0 only if y Q = 0. Error returns are: 

IERRF = 0; normal return 

IERRF = 1; the leading coefficient is zero 

IERRF = 2; one of the roots failed to converge in the initial iteration cycle 

IERRF = 3; one of the roots failed to converge in the improvement iteration cycle. 

INTI A is a closed subroutine for the solution of a set of simultaneous differential 
equations. It is a variable- interval- size routine and uses the fourth-order Runge-Kutta 
formula in conjunction with Richardson's extrapolation to the limit theory. Error returns 
are: 

IERR = 1; normal return 

IERR = 2; ELT block is not monotonic in the direction of integration 
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IERR = 3; variables have failed to meet the local truncation error requirements 
nine consecutive times 

IERR = 4; variables have failed to meet the local truncation error requirements 
at least nine times over the last three intervals. An acceptable 
answer has been reached, however, and is in the VAR array. 

ITR2 finds a value for x within a given epsilon of relative error in a given interval 
for a given F(x) = 0. Error returns are: 

ICODE = 0; normal return 

ICODE = 1; maximum iterations (= 150) are exceeded 

ICODE = 2; DELTX (the scanning interval) = 0, or negative 

ICODE = 3; a root cannot be found within the given bounds, ALLU and BULU 

ICODE = 4; ALLU > BULU 

CHSUB and DERSUB are subroutines required by INTI A 

FOFX is a function called by ITR2. 


Options 

An option is available for including or omitting (3 from the X equation. A value 
for /3 can be read in and will be included in equations (23), or it need not be read if it 
is not to be included in the X equation. The variable controlling this option is listed 
under "Input." 


Input 

The NAMELIST statement is used to put data in. NAMELIST names and the variable 
names contained in each with their explanations are listed in the following table: 


NAMELIST 

name 

Variable names 

Explanations of variables 

NAM1 

A 

A 


SIGMA 

a 


ALPHA 

a 


PRMIN 

The absolute value of an increment of the indepen- 
dent variable which is the frequency of printing 
results. (See "Output.") 
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NAMELIST 

name 


Variable names 


CH, CIMAX 


Explanations of variables 

The initial computing interval and maximum desired 
computing interval, respectively, as required by 


ALLUI, BULUI, 
DELTX 


EP1, EP2 


INT1A. 

The lower bound or initial guess for u, the upper 
bound or final guess for u, and the size of the 
scanning interval, respectively, as required by 
ITR2 in solving for u. 

Relative and absolute error criterion, respectively, 
used in ITR2 to stop the iteration when either of 
these convergence criteria are satisfied: 


1. If | Ui |>EPl, 


u i ~ u i-l 
u i 


^ EP1 


NAM2 


2. If u t ^ EP1, 

| u i " u i-! I = EP2 


NT 

IOPBETA 


VARO 


Values of EP1 = EP2 = 0.1 x 10 -6 have given 
satisfactory results. 

The number of values in the ELT block described 
below. 

Equals 0, j3 is not read and is not included in 
the X equation. 

Equals 2, (3 is read and is included in the 
X equation. 

A one-dimensional array containing the initial 
values of the independent variable followed by 
the three dependent variables, x, y, and X, 
in that order. 
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NAMELIST 

name 

Variable names 

Explanations of variables 


ELEl, ELE2 

One-dimensional arrays containing the upper bounds 
of relative truncation error and "relative zeros," 
respectively, for the dependent variables as 
required by INT1A. If the error for any variable 
exceeds its respective ELEl value, the computing 
interval is halved and the integration restarted at 
the beginning of the interval. Under certain cri- 
teria the interval is doubled. If the absolute value 
of any of the variables is less than its respective 
ELE2 value, the relative error criteria for that 
variable will not be applied. Satisfactory results 
have been obtained using ELEl values = 10"^ 
and ELE2 values = 10 - ®. 


ELT 

One- dimensional array of NT values, monotonic in 
the direction of integration, at which the user 
specifically desires control returned to his pro- 
gram from INTI A. (See "Output.”) 

NAM 3 

BETA 

i 

( 3 , read only if IOPBETA = 2. 


Output 

The frequency of printing results is determined by the input quantity PRMIN. The 
initial values are printed and the results will be printed again when the independent value 
is first greater than PRMIN and thereafter when the independent value has been updated 
by at least as much as PRMIN. The final value in the ELT array should be the final value 
of the independent variable (equals 1.0). INT1A iterates to obtain results at the specific 
values listed in ELT. If results at values of x other than 1.0 are desired, the program 
can be easily modified to have them printed. Values of the drag coefficients are printed 
at the end of each case. The total drag coefficient is given, as well as the pressure and 
friction drag coefficients. 


Computational Flow Diagram 

A concise computational flow diagram is included here to show the steps in com- 
puting. Details can be readily obtained from the listing. 
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Listing of FORTRAN Program 

The FORTRAN program, including FALG, INT1A, and ITR2 (Langley Research 
Center library subroutines) is as follows: 

PROGRAM BODREV ( INPUT# OUTPUT , TAPB5= I NPlJT , TAPF6 = OUTPUT ) 

C PROGRAM TO COMPUTE BODIES OF' REVOLUTION HAVING MINIMUM VISCOUS 

C PLUS (NEWTONIAN) PRESSURE DRAG IN HY^PSOMIC FLOW* 

C VARIABLES NEEDED IN INTEGRATION ARE STORED IN THF VAR ARRAY 

C (TENTATIVF ANSWERS IN CORRESPONDING POSITIONS IN THE CUVAR 

C ARRAY) 

C VAR ( 1 ) INDEPENDENT VARIABLE 

C VAR (2) X 

C V A R ( 3 ) Y 

C VAR (4) LAMBDA 

C 

C Y PRIME IS REFERRED TO AS U IN THIS PROGRAM 

C 

C OPTIONS 

IF IOPBFTA = O, B C T A IS NOT RF AO AND IS NOT INCLUOFD IN THF 
LAMBDA DOT FOUATION* 

IF IOPBETA = 2, BETA IS READ AND IS INCLUDED IN TH^ LAMBDA DOT 
EQUAT I ON* 

COMPLEX ROOTS*TF(vp 

D I MENS ION ELE1 < 3 > *ELF2 < 3 > « ELT ( 1 0 ) « ERRvAL ( 3 ) * COEFFS (1 0 ) * ROOTS ( 4 ) , 

1 TFMP (10) 

EXTERNAL DERSUH fCHSUB 

COMMON /BLK 1 /CUVAR ( 4 ) , VAR ( 4 ) , OER ( 4 ) * VARO (4 ) ,A, SIGMA, ALPHA * U * ALLIJ * 
1 BULU* DELTX * fpi ,FP2* I C OOF * I I , C I , BFT A , I 0 D BFT A , I OPUO « l.JO 
NAMFL I ST /NAM! / A# S I CMA, ALPHA , PRM I N * C I I * C I MAX * ALL U I , B'JLU I * DELTX * 

1 FP1,FP?,NT*I OPRFT A 

2 /NAM2/VAR0*FLF1 *FLF?, C LT 

3 /NAM 3/RET A 

0 READ (SiNAMl ) 

READ ( E , NAM? ) 

READ QUANTITIES IN NAM3 ONLY IF THFY ARF REQUIRED FOR THE. RIJN. 

SEE I OPUO AND IOPBETA AROVF • 

READ (5 » NAM3 ) 

NCASE=0 

IF(VAR0(3)*FQ*°.0)NCASF=1 

1 C I =C I I 
PFR=^#999R9QQQ 
N= 3 

I TEXT = 0 
SPEC = 0 % 0 

I I 

nplusi =N+1 
DO 20 1=1. NPLUS 1 

0 VAR ( I ) = VARO ( I ) 

allu= ALLU I 
BULU = BULU I 
WR ITF ( 6 ♦ NAM 1 ) 

WRITE (6,NAM2) 

WRITF (6*NAM3) 

WRITF (6,70) 

WRITE (6,80) VAR ( 3 ) , A , S IGMA , ALPHA 
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I F ( VARO (3 ) #NF .0 .0 ) WR I TF (6*90) 

IF ( VARO (3 ) • F Q • O . O ) W R I T F (6i 100 ) 

IF (IOPBFTA.FQ.O) WRITF (6*120) 

IF ( IOPBFTA.BO.?) WRIT 11, (6*130) BETA 
C 

C IF Y IF INITIALLY = 0, THE INITIAL U WILL BE OBTAINED FROM THE 

C QUART I C IN U <I0PU0 = 2>* THERE WILL BE TWO REAL ROOTS* UOl 

C AMD U02. COMPUTATIONS WILL BE DONE FOP ROTH CASES. IF Y IS 

C NOT = 0 INITIALLY* THE INITIAL U IS OUT A I NED FROM THE 

C IMPLICIT EULER EQUATION ( I OPUO = 0)* 

I F ( V ARO (3) .NE.O.O) I OPUO =0 
IF ( VARO (3 ) .N".0.0 >G0 TO RE 
IF (NCASE.NE.2 )GO TO 26 
U0=U02 
GO TO 22 
26 CONTINUF 

I OPUO = 2 
NDEG=4 
I C OFFF = 0 
COEFFS (1 ) =1 .0 

COFFFS ( ? ) =-4 .0* ( S I GMA*#ALPHA )/A 
COEFFS (3 > =2. O 
COFFFS (4 )=0.o 
COFFFS (F) = l .0 

CALL FALG (CO-^FS * ND~G * I COE FF , ROOT^ * TFMP* IFRRF ) 

IF ( IERRF.EO.OGO TO 23 
WR I TE ( 6 * 1 60 ) 1 ERPF 
GO TO 10 

23 CONTINUF 
I RR= 1 

no ?4 1=1*4 

R I =A I MAG (ROOTS ( I ) ) 

IF (RI .NE.O.O )GO TO 2* 

IF ( IRR.EQ.2 )G0 TO 25 
UO 1 =REAL (ROOTS ( I ) ) 

IRR = 2 
GO TO 24 

25 U02 = REA L (ROOTS ( I ) ) 

24 CONTINUF 

uo=uo 1 

2 CONTINUF 

initialize in subroutine intia 

CALL INTIA ( I I ,M,NT*CI * SPEC *C1 MAX* I ERR * V AR * C U VAR * D^R . EL F 1 *FLE2*FLT 
1 , srrval*deRsur*chc U r* ITFXT ) 

C COMPUTATION FOR DRAG INTEGRAL (INITIALIZATION) 

C DRAG 1 CONTAINS THE Y0**2 TERM PLUS THE FIRST PART OF THE INTEGRAL 

C DRAG? CONTAINS ONLY THE A TERM OF THE INTEGRAL 

FD1 M 1 = VAR ( 3 ) * ( 2 . 0*U**3/ ( 1 • 0 + lJ**2 ) ) 

FD2M1=VAR(3 ) * ( A/ ( (SI GMA+VAR ( ? ) *SGRT ( 1 • 0+U**2 ) ) **ALPHA ) ) 

DRAG 1 = V AR ( 3 ) *-*2 
DRAG 2=0.0 
DR AG = DR AG 1 +DRAG2 
ST ORFT = VAR ( 1 ) 

WRITE (6*140) 

GO TO FO 

C start integration 

30 CALL INTIA (I I * N * NT * C I *SPEC*CIMAX* I ERR * VAR * CU VAR , DER ♦ ELF 1 * ELE2 * ELT 

1 *FRRVAL*DERSUB*C.HSUB* ITFXT ) 

IF ( ( IERR.EQ. 1 ) .OR. ( IERR.EQ.4 ) ) GO TO 40 
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WRITE (6*110) I 
GO TO 10 
40 CONTINUE 

IF (II •EQ.2 ) GO TO 30 
C COMPUTATION FOR DRAG INTEGRAL 

FD 1 =VAR ( 3 )* ( 2*0*U**3/ ( 1 .0+U**? ) ) 

FD2=VAR(3)*(A/( (SIG^A+VAR(2) *^QRT ( 1 • 0+U**2 ) ) ^^ALPHA ) ) 

FDA i = (FD1 Ml + FD1 )/2.0 
FDA2= (FD2M1 T.F02 >/2.0 
FD I NT 1 = FDA 1 * ( VAQ ( 1 ) — Q T ORET ) 

FDINT2=FDA2* ( VAR ( 1 l-STORET ) 

DRAG 1 = DRAG 1 +FD I NT 1 
DRAG2 = DRAG2 + cr Q I NT 2 
DR AG=DRAG 1 +DRAG2 
FD1M1 =FDi 
FD2M 1 =FD2 
STORFT = VAR ( 1 ) 

C TEST TO SEE I f VALUES SHOULD RF PRINTED 

IF (VAR( 1 ) •GF.Pcp ) GO TO SO 
PRFRFQrrVAR ( 1 ) — PPV AL 

IF (PRFREQ.GTtPpMlN) GO TO RO 
IF ( VAR ( 1 ) tLT.P^R ) GO TO 30 
SO CONTIMtF 

PRVAl=VAR ( I ) 

WRITE (6*60) VAR (2 ) * VAR (3 ) * U 
IF ( VAR ( 1 ) .LT.PFR ) GO TO 30 
C0NVD=2*C/(VAR(^)**2) 

CDT=DRAG*CONVD 

CDP- DR AG 1 *C ONVD 

CDF=DRAG2*C0NVD 

WRITF (6*150) COT *CD D *CDF 

IF ( IOPUO.FQ. 0)00 TO 10 

IF (NCASF.EQ*? ) GO TO IO 

NCASF=P 

GO TO 1 1 

c 

60 FORMAT (4F18.8) 

70 FORMAT ( 1 H 1 1 X4 ^HROD I f S OF REVOLUTION HAVING MINIMUM V I S COU S 1 X44 H PL 

1US (NFWTONIAN) PRESHjQ^ DRAG IN H YPER^ON I C 1 X4 H^LOW ) 

80 FORMAT ( /// 1 PX5H INPUT 16X2HYQ 1 7X 1 H A 1 3 XSHS I GMA 1 3X C >H AL P H A /E 36 • 8 * 3E 1 8. 

1 8 // ) 

on FORMAT ( 1 8 X* INITIAL U OFT A I N^O FROM FtjLER FGUATION#) 

100 FORMAT ( 1 8X# INITIAL J OBTAIN^O FROM QljARTIC EQUATION*) 

1 IO FORMAT (1X32HFRDOR RFTljRN FROM INTI* I ERR = 14//) 

120 FORMAT (18X40HPFTA NO T INCLUDED IN LAMBDA DOT EQUATION) 

130 FORMAT ( 1 8 X46HBFT A INCLUDED IN LAMBDA DOT EQUATION. BF T A =E16*8) 

140 FORMAT (//14X1HX17X1HY1 1X7HY PRIME/) 

ISO FORMAT (//6X1 6HT0TAL DRAG C0FFF7X 1 5HPRF S DRAG COFFF 3 X 1 9HFR I C T I ON D 
1 RAG C0FFF/3F22.8 ) 

160 FORMAT ( 1 X*ERROR RETURN FROM F ALG • 1 FRRF = *-,14) 

END 
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SUBROUT I NF O r RF| jR 
EXTFRNAL FOFX 

COMMON /BLK 1 /CUVAR ( 4 ) * VAR < 4 ) • HER (4 > ,VAR0(4) f A 1 SIGMA 1 ALPHA * U * ALLU « 

1 P ULU n DEL TX * 1 I CODF * I I , C I * BET A ♦ I OPPFTA , I OPUO ♦ UO 

IF ( ( IOPUO.NF.p ) .OR* ( I I .NE.O ) ) GO TO 1 0 
U=UO 

GO TO 30 

10 CALL I TR2 (U * ALLU * BULU * DEL TX ♦ F OFX i £P 1 *EP2n 1 50 * I CODF ^ 

IF (ICODF.EQ.nj GO TO 20 

FUNU = CUVAR (3 >* ( ?»0*U**2* < 3 • 0+U**2 ) / ( (1 • 0+U**2 ) ) - A LPH A *A *CU V A R (2 
1 )#U/( (SORT ( 1 • n + l J**2 ) ) # ( (SI GMA+CUVAR ( 2 )*SGRT (1 *0+U**2 ) )** ( 1 .0 + ALPHA 
2 ) ) ) ) +CUVAR ( 4 ) 

IF (FUNU.LT.n. ic-06) OO TO 20 
WRITE (6 #40) ICODFtU*FUNU 
STOP 

20 CONTINUE 

3n C0N7 INUF 

DFR (1 ) =0 . 0 
DER ( 2 ) = 1 .0 
DFR (3 ) =U 

DFR (4 ) =-?.0*U**3/ ( 1 ♦ 0+U**2 ) -A/ ( ( S I GM A+CU VAR ( 2 >*SORT ( 1 • 0 +U* *2 ) ) ** AL 
1 PHA ) 

IF (i OPPFTA • "O • P ) n~R ( 4 ) = D cr R ( 4 ) -2 • 0 #RfTA*CIJVAR ( 3 ) 

RETURN 

r 

40 FORMAT ( 1 X 3 3 HE POOR RETURN P ROW I TP2 . ICODE = I4,3x4HU = F16*8t3X7 
1HF(U) = F i 6 . B// ) 

END 


FUMCT I ON F OF X ( X ) 

COMMON /RL«1 /CUVAR ( 4 > «.VAR < 4 ) ♦ DFR ( 4 ) , V A R O ( 4 > , A , S I GM A , ALPHA « U « ALLU ' 

1 BUL U » DELTX, fpi , EP2 ♦ I COD- ,1 HCHB^TAH OPRETA , IOPUO ♦ UO 
F0Fx = 2 • 0*X**2* ( 3. 0+X**2 ) / ( ( I • D + X**2 ) **2 ) - ( ALPHA*A*CUVAR < 2 ) *X ) / < ( SQ 
1 RT ( 1 #0+X**2 ) )* ( ( S I GMA + CUVAR ( 2 )*SGRT ( 1 • 0 + X**2 ) ) ** ( 1 • 0 + ALPHA ) ) ) +CUVA 
2R(4 ) /CUVAR (3 ) 

RETURN 

FMD 
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SUBROUT INE F ALG ( COEFFS , N i I *ROOT «TFMp,IERR) 

-H-****#*** nOCUMPNT hajp 08-01-68 SUBROUTINE RFVISED 08-0 1-68 ********* 
DIMENSION COEFFS (1 ) » TF MP ( 1 ) .ROOT ( 1 ) 

D I MENS ION X I D I FF < 2 > • RD I FF < 2 > , APPROX ( 3 ) 

COMPLEX F.FPR, APPROX, TEMP .ROOT 

complex tempm 
COMPLEX RELTST 
NSA VF=N 
I ERR = 0 
I B = 0 

I CLEAN=2*N+2 

C CLEAR WORKING Area 

f DO 7771 LLL =1*1 CLEAN 

7771 TEMP (LLL ) =0* 0 

c clear root storage 

DO 7772 LLL = 1 , N 
> 7772 ROOT (LLL )=0.0 

CONSTANTS to test 
CONVERGENCE 
CONST “ • 1 E— 6 
C OVERFLOW 

OVCON= 1 .F150 

C MAGNITUDE OF ROOTS 

RCONST = 1 .F-21 

C JONJON = 0 • F I RST ITERATION 

JON JON = 0 

C CHECK CONSTANT TFRM FOR 7FRO 

JJJ=1 
NCO=N+l 

802 IFd.NE.l )GO TO 800 
C COMPLEX COEFFICIENTS 

NC0=2*NC0 

IF (COEFFS (NCO-1 ) .NE.O. ) GO TO 101 
C HERE IF R^AL COEFFICIENTS 

800 I F ( COEFFS ( NCO ) • N r • 0 • ) GO TO 101 

C ROOT=ZERo 

801 ROOT ( J J J ) = 0 • 

NC O ~ NCO— 1 
JJJ=JJJ+ 1 

C REDUCE DEGRFF AND IF 1 » ST ORE ROOT AND EXIT 

N = N— 1 

IF (N.NE. 1 ) GO TO 80? 

ROOT ( J J J ) = 0 • 

GO TO 1 006 
C 

C ENTRY FIRST AND SECOND ITERATIONS 

PI J = JJJ 

NT ERMS = N+ 1 
KC ON J = 0 

C CLEAR APPROX 

APPROX (1 ) =0*0 
APPROX ( 2 ) =0 • 0 
APPROX ( 3 ) = 0 • 0 
IE ( I • EQ • 1 ) GO TO 43 

C REAL COEFFICIENTS 

DO 78 I JFF= 1 .NTFRMS 

78 TEMP ( I JFF )=CMPLX (COEFFS ( I JFF ) ,0.0) 

GO TO 700 

C COMPLEX COEFFICIENTS 

43 DO 79 I I IX=1 . NTFRMS 

79 TEMP ( I I IX ) =CMPLX (COEFFS (2*1 I I X-l > * COEFFS (2* I I I X ) ) 
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C CHECK LEADING COEFFICIENT FOR O OR 1 

700 TEMPL=REAL (TEMP ( 1 ) ) 

I F (TEMPL.NE. O. >G0 TO 701 
C IF REAL IS ZERO, CHECK IMAGINARY 

TEMPL = A I MAG ( TEMP ( 1 ) ) 

IFCTEMPL.NE.O* )G0 TO 701 

C LEADIN'^ r O tr ' rfr I r T C NT 7 C RO I^RR=1 

I F P P — 1 
GO TO 1006 

C DIVIDE BY LEADING COEFFICIENT 

701 TFMPM=TEMP ( 1 ) 

DO 70? LL A = 1 fNT^PMF f 

70 p TFMP (LLA ) =TFmd ( L | A)/TFMPM 
C KC ON J = 1 i TRIAL VALUE= CONJUGATE 

c 

746 IF (KCONJ.NE. 0 )G0 TO 47 
C FIRST TRIAL VALU^ 

A D PROX ( 1)= ( *01 , .01 ) 

C DIFFERENTIATE 

47 DO 8 11=1 ,NTFRMC 

xpon=nterms~ I I 

NNOW = I I+NTERMS 
8 TEMP (NNOW ) =XPON*TFMP (II) 

NPON = NTERM S — 1 

C K A = O FOR FIRST TRIAL VALUE 

K A = O 

C 

C JON JON = 1 SECOND ITERATION 

749 IF ( JONJON.EQ. 1 )A°PR0X(1 ) =ROOT ( J ) 

c clear rhiff.xi'Mff 

DO 7773 LLL =1*2 
RDIFF(LLL) “0*0 
7773 XIDIFF(LLL )=0.O 
C 

C ROOT EVALUATION 

C MAXIMUM ITERATIONS =120 

13 L = 2 

PARTRl = REAL (APPROX ( 1 ) ) 

P ARTM 1 = A I MAG ( APPROX ( 1 ) ) 

DO 1 ? K=? , 1 2 ) 

C EVALUATE F(X). 

10 E= (0.0, o,o ) 

DO 9 11=1 ,NTEQMS 

F = APPR0X(L-1 >*F+TEMP( I I ) 

XF = ABS (REAL (E ) ) 

YF = ABS ( A I MAG (F ) ) 

C CHECK FOR OVERFLOW 

I F ( XF •GT.OVCON.OR. YF *GT*OVCON)GO TO 14 
9 CONTINUE 

C EVALUATE FPRIME(X) 

FPR= (0.0, 0.0 ) 

DO 1 1 JJ=1 , NPON 
NNOW= JJ+NTERMS 

FPR= APPROX (L-l ) #FPR+TEMP ( NNOW ) 

YFP-ABS ( A I MAG (FDR ) ) 

XFP=ABS (REAL (FPR ) ) 

C CHECK ^OR OVERFLOW 
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l 1 

c 

c 


c 


c 

6732 

C 

732 


18 


87^ 


1 2 


r 

1 4 


C 

1 36 
C 


I F (XFP. GT • OVCON. OR • YF^.GT.OVCON )GO TO 1 4 
CONTINUE" 

SEE IF FPRIME=0 

IF (XFP.EQ.O. 0. AND • YFP • C ’Q*0*0 ) GO TO 14 
IF NOT ZERO « NEW APPROXIMATION 
APPROX ( L ) = APPROX ( L— 1 ) -F/FP R 
PA RTR2 = REAL ( APPROX ( L ) ) 

PARTM2 = AIMAG ( APpROX ( L ) ) 

SFT F I T HER PART TO ZERO IF LFSS THAN l.F-21 
IF ( A8S ( PARTR2 ) .LE.RCONST )PARTQ2=0. 

IF (ABS (PARTM2 ) .LE.RCONST )PARTM2=0. 

I F ( PARTR2 • EQ • o • •AND.PARTM2.FQ, 0* ) GO TO 6732 
GO TO 732 

ZERO ROOT 

IF (L .EG# 3 ) APPROX (2 ) = APPROX ( 3 ) 

GO TO 81 

RD IFF ( L— 1 ) =ABS (PARTR2-PARTR1 ) 

XIDIFFIL-1 ) = ABS ( PARTM2— PARTM 1 ) 

IF(L.EQ.3) GO TO 18 
L = 3 

PARTR1 =PARTR2 
PARTM 1 =PARTM2 
GO TO 1 0 

TEST 1 

I F ( (RD I FF (2 ) +X 1 O I FF (2 ) ) *LT • (RD I FF ( 1 ) + x I D IFF ( 1 ) ) ) GO TO 
TEST 2 

RELTST = (APPROX ( 3 ) -APPROX ( 2 ) ) / APPROX ( 3 ) 

D I FFR = A6S (REAL (RELTST ) ) 

D I FFX I =ABS ( A I MAG ( RFLTST ) ) 

I F ( D I FFR.LT • CON FT. ANO.D I FFX I .LT .CONST ) GO TO 8 1 
APPROX ( 2 ) =CMPLX (PARTR2 , PART M2 ) 

PARTR1 = PA RT R 2 
PARTM 1 = PA RT M 2 
R D I FF ( 1 )=RDIFF (2 ) 

X I D I FF ( I ) = X I O I ff* ( 2 ) 

MAXIMUM ITERATIONS EXCFFDFO OR 
OVERFLOW OR 
FPR I ME = 0 

TRY AGAIN WITH SECOND TRIAL VALIF 
I F ( JONJON.EQ • 1 ) ^O TO 136 
IF (K-A.erO, 1 05 )00 TO 7l 
APPROX { 1 > = ( 1 • , 1 * ) 

K A = 1 05 
GO TO 1 3 

SECOND ITERATION NONC ONVERGENT ROOT IFRR=3 
I ERR=3 

STORE RFSULT AND I MPROVF NEXT ROOT 
GO TO 82 


8700 
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FIRST ITERATION ROOT R NONCON VERGF NT IFRR-2 
IMPROVF (R-l ) ROOTS 
71 IFRR=2 
C I 8 =L AST CONVERGENT ROOT 

I B = J 

ROOT ( J ) = APPROX ( P ) 

IF ( IB*NE. 1 ) GO TO 9971 
C FIRST ROOT p A I LED RETURN 

GO TO 100* 

9971 JON JON = 1 
GO TO 101 
r 

C STORE ROOTS 

fll I e ( JON JON • EQ • 1 ) GO TO 82 

C HERE IE FIpST ITFPATION 

5 ROOT ( J )=APPROX (P ) 

C RFDUCE POLYNOMIAL RY SYNTHETIC DIVISION 

1 NTERMS=NTERMS- 1 
DO 7 IK=2,NTERMC 

TEMP (IK) =ROOT ( J > ^ TEMP ( I K— 1 )+TEMP (IK) 

7 CONTINUE 

C NEXT ROOT IF COMPLEX COEFFICIENTS 

IF ( I .EG. 1 ) GO TO 745 
C HERE IF P^AL COEFFICIENTS 

I F ( KCONJ • FQ* 0 ) GO TO 744 

r RESET KCONJ IF POOT IS CONJUGATE OF PREVIOUS ROOT 

KC0MJ=0 
GO TO 745 

744 X=RFAL(ROOT (J) ) 

Y= A I MAG (ROOT ( J ) ) 

IF ( X • EQ • 0 • ) GO TO 745 
C SEE IF REAL OR COMPLEX 

I F ( ABS ( Y/X ) *LE • 1 . E- 1 0 >G0 TO 745 
C COMPLEX root trial valuf=conjugate 

APPROX ( 1 ) =CONJG (ROOT ( J ) ) 

KCONJ= 1 

C NEXT ROOT 

745 J = J + 1 

IF ( J.NE.NSAVF )GO TO 746 
C LAST ROOT 

ROOT ( J ) =-TEM^ ( 8 ) /TEMP ( 1 ) 

JON JON= 1 
GO TO lOl 

C 

C IMPROVFD ROOT FROM SECOND ITFRATION 

es ROOT ( J )= APPROX (2 ) 

X=ABS(REAL (ROOT ( J ) ) ) 

Y = ABS ( A I MAG (ROOT ( J ) ) ) 

C SET REAL OP IMAG. TO ZERO IF LESS THAN RCON^T 

IF ( X • LT • RCONFT • ANO. y # LT • RCONST ) GO TO 53 
IF ( X • LT • RCONST ) ROOT ( J ) =C^LX (C.O*AIMAG(ROOT (J) ) ) 

I F ( Y*LT .RCONST ) ROOT ( J ) =CMPLX (RFAL ( ROOT ( J ) ) • n. n > 

GO TO 108 

33 IF(X.G r *Y) GO TO 34 

ROOT ( J ) =CMPLX ( 0 • 0 # A I MAG ( ROOT ( J ) ) ) 

GO TO 108 

34 ROOT ( J ) =CMPLX (PEAL (ROOT ( J ) ) i 0 * 0 ) 
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c 

1 n 8 
C 

1 006 


1 12 
1 14 

1 

1 004 

l oon 

4 

999 

6 

7 
3 

1 003 

1 002 

5 

998 

9 

8 


IF(R-l) ROOTS I MPROVEO , RETURN 
IFfJ.EQ. IBJGO TO 10C5 
J= J+l 

IF N ROOTS I WDROVFD * RETURN 
IF(J.LF*NSAVF)GO TO 749 
N=NS AVE 
return 
end 


SUBROUT I NE I TR2 ( X , A , R , DFLTX,FOFX» El »F2 iMAXI , I COOF ) 
Y= A 
KX=0 
LX = P 

I F (DELTX >111,111,112 
IF <B- A ) 113,113*114 

I =0 

IF ( F OFX ( A ) ) 1 ,2,3 
XB 1 = X 

IF (LX.N'F.O >GO TO 1 00 1 
X =X+DFLTX 

I F ( X — B >1000, 1000,1004 
X = B 
LX = 1 

I F (FOFX ( X ) ) 1 , 2 ,4 
XB = X 

X = X— DEL TX/ (2*** ( I +1 ) ) 

1 = 1 + 1 

IF (MAXI ,LT* I )G0 TO 444 

I F (FOFX ( X ) >6,2,7 

L = 1 

XX-XR 

GO TO 18 

L = 2 

X X = XB 1 

GO TO 1 8 

XB 1 = X 

IF (KX*NE,0 >G0 TO 1001 
X = X + DELTX 

I F ( X-B >1002,1002,1003 
X = B 
KX= 1 

IF(FOFX(X) >5,2,3 
XB = X 

X=X-OELTX/ (2.** ( I +1 > > 

1 = 1+1 

IF (MAXI .LT. I >GO TO 444 

I F (FOFX ( X ) ) 8 2 , 9 

L = 3 

XX = XB 

GO TO 18 

L =4 

XX = XB 1 
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18 IF ( A3S (X ) —E 1)36* 36 *37 
37 IF ( ABS < ( XX— X ) /X > — F 1 ) 2 * 2 * 1 7 
76 IF ( A5S (XX-X )-F? )2*2* 17 
17 GO TO (81 « 4 « 8 1 * 5 ) * L 
81 X81 = X 

X=X+nF|_TX/ (2«** ( I +1 ) ) 

GO TO (999*4*998*^) *L 

1 1 1 ICODF = ? 

GO TO 79 
1 1 3 ICODF =4 
GO TO 79 
r°] ICODF =3 
GO TO 79 
444 I CODE = 1 

GO TO 79 
2 ICODF =0 
79 CONTINUF 
RETURN 
FND 


SU6R0UT I NE I NT 1 A (I I *N*NT *CI ♦ SPEC * Cl MAX, I ERR * V AR , CUV AR * DER ♦ ELE 1 
1 FLEE , EL T * ERR VAL. * DFRSUB , CHSUB , I TEXT ) 

D I mens ION SI VAP (20)* SELF 1 (20 > * FLF 1(20)* ELF 2 ( 20 ) * DEP (21)* 

I FOFPv (21 ) * ROY ( 20 ) * SDY 1 ( 2 0 ) , Y I MCP (20 ) *FPPVAL (20 ) * ERVOVH (20 ) 

? C LT (10), S r L T (13)* R^LM T N ( ?0 ) * STFP ( 3 ) 

D I MFNS I ON VAR ( ?l > * C U V A P (21 ) 

I NTFGER TEX (15) 

I NTEGER CODE , TP^H , SUMHAF ♦ STEP , TEST , DCODE 
RFAL K 1 

BEG IN IN I T I AL I ZAT I ON 
IF (II *GT. 0 > GO TO 520 

1 TP = T 

S S P E C = S IGN ( SPEC *C I ) 

SCIMAX = SIGN(CIMAX*CI ) 

VAR1 = VAP ( 1 ) 

IF (Cl • F O • 0*0) GO TO 151 

IF (SSPEC • f Q # o.O) GO TO 7 

IF (ABS(SCIMAX) • GT • ABSCSSPFC ) • OR • SCIMAX #EG§. 0*0 ) 

1 SC I MAX = SSPEC 

TEST TO SEE IF vAR IS ZERO 

IF ( ABS ( VAR1 ) *GT . 1 .OE- 1 1 )G0 TO 2 

TP=SSPFC 

GO TO 7 

2 IF ( (VAR1 /SSPEC ) .GT. 1*F>13) GO TO 4 

3 K1 =C.O 
G0T06 

4 <1=1*0 

6 TP = VAR 1 -A MOD ( VARl , SSPFC ) 

I F ( ABS ( r°-VARl ) *LT. 1 • F- 12X1=1.0 
TP=TP+K1*SSPEC 

IF ( ABS ( (TP-VARl )/VARl ) • L T • 1*E-11> TP = TP + SPEC 

TEST FOR DIRECTION OF INTEGRATION 

7 <1=1,0 

IF (Cl • L f • 0*0) <1=-1*0 
CIK=CI*K1 
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C I MAxK=SC IMAX*K1 
TPK = T P*K 1 
\/ARK = VAR 1 *< t 

C SET UP STORAGE POP INTERNAL U Q E 

NP 1 = N+1 
NELT = 1 
REMA I N=0 • 0 
MHAF= n 
MTS=mT 
SUMHAF = 0 
LOOP = O 
DO 91 1=1 *3 

91 STEP ( I ) =0 
I ERR= 1 

DO 8 1=1 » NP 1 

8 CUVAR ( I ) = VAR ( I ) 

DO 101 1 = 1 iN 

1 01 SELF1 ( I ) =ELF_'l ( I ) 

IF (NT • EO • O ) GO TO 13 
ICO IF (NT .EG* 1) GO TC 10 
NT M 1 =NT— 1 
FLTK=K1 *flt ( 1 ) 

DO 9 1 = 1 ,NTM 1 

ELTK2=K1*FLT ( 1 + 1 > 

IF (FLTK • LT • FLTK?) GO T09 
GO TO POO 

9 FLTK = EL TK 2 
1 n CONTINUE 

ELTK=K1 *F|_T ( NFL T ) 

IF (VARK • L T • <=LT<) GO TO 1 1 

IF (NELT .EG. NT) GO TO 1 3 
NELT =NEL T+ 1 
GO TO 10 

1 1 NFLTL=NT-NPLT+1 
GO T012 

13 NEL TL = 0 

12 DO 14 I = 1 « N 

14 RELMIN( I )=SFLE1 ( I )/l 28*0 

IF (NT .PQ. 0) GO TO 906 
DO 995 1=1 (NT 

995 SELT ( I )=ELT ( I ) 

996 CALL DERSU3 

IF (II • F O • 4 ) GO TO 120 

DO 15 I = 1 f N 

1 5 FDERV ( I ) =DER ( I + 1 ) 

11=1 

TEST =0 

DO 300 1 = 1 * 1 P 

300 TFX ( I ) =0 
TEX ( 1 ) = 1 
TEX (2 ) = 1 
KK 3= 1 

IF (ITEXT) 635 % 63 ♦ 6 35 
151 PRINT 1000 

1000 FORMAT (//11H Cl IS 7FR0) 

STOP 

C END OF INITIALIZATION 
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520 11=1 

T P S H = 0 
LTSH=0 

VARK=VAR (1 >*K1 
C I K=C I *K 1 

S 1 =vark+c jk 

IF (SSPFC • FO* O.O) GO TO 52F 
Kl<= 1 

IF (NELTL.EQ. 0) GO TO 17 
IF ' (FLTK-TP K ) 16*16*17 

16 CV = P(_ TK 
COOF= 1 
GO TO 1 8 
]~7 CV-TPK 
^ n n cr = p 

18 IF (ABS CCV) .LT. 1 # c -12 ) GO TO 680 
I f f CV -fj )20 ,p^ , i q 

lO I F ( ARC ( ( CV~F 1 ) /r\i ) , . . 1 F- 1 1 ) GO TO 555 

pn if (NFLTL • er Om O) GO TO 540 

IF(ARS( (FLT<-T^)/CV).LT..1?-1 1 )GO TO 550 
IF (CODE .EQ. l) GO TO 545 
54^ DX ~ TP — VAR ( 1 ) 

T F X ( 5 ) = 1 
TP=tp+ SSPEC 
TPK=TP*K 1 
TPSH= 1 
GO TO 560 

C FHORT INTERVAL nu- TO BOTH 

5 cr '' TP = TP + 5FPFC 
TFX (6 ) = 1 
TPK=TP*K 1 
f PSH- i 
GO TO 545 

C IF HERE CV IS LIKELY ZFRO 

530 IF (SI .LT.-l .OE-12 ) GO TO 536 
IF (CODE .EQ, 1) GO TO 550 
IF (NELTL • E O • O) GO TO 540 
IF ( ABS (FL TK ) .LT. 1 •OE-12 )G0 TO 550 
GO TO 540 

C SPEC IF 7^P0 

F P f I F ( ABS ( RFjvA I N ) .OT • • 1 cr- ] 1 ) GO TO 06 
IF (NFLTL . c 0. o) GO TO 565 
94 I F ( ARS ( ELTK ) .GF . 1 .F- 1 2 ) GO TO PI 
IF (SI .LT.-l. OF- 1 2 ) GO TO 566 
GO TO 545 

21 S2=FL1K-S1 

•IF (S2 ) 545*545*22 

22 I F ( A8S ( S2/ELTK ) .LT. 1 . QF-1 2 ) GO TO 545 
GO TO 565 

C shoptintfrval If DUE TO ELT block 

^45 Df L T= SFlT(M-LT) 

TFX (4 ) = 1 

DX = DFLT- VAP ( 1 ) 

RFMA IN=C I -DX 
RFMA I K=RFMA I N*F 1 
LTSH= 1 

nelt=nflt+ 1 
NFLTL = NE L TL - 1 
IF (NELTL.EQ.O)GO TO 560 
FL TK -K 1 *SFLT (^| T ) 

GO TO 560 
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565 0 X = C I 

TFX (3 ) =1 
GO TO 560 

96 IF (NFLTL *EQ, O ) GO TO 98 

IF (FLTK *LT • < \/AR«4 RFm A I < ) ) GO TO 9* 

98 DX = RFM A I N 
TEX (7 ) = 1 

rfma i n= o • o 

C-O TO 560 
535 DX = C I 

TEX ( 3 ) = 1 
TEST = 1 
GO TO 555 

QFG I N RUNGE-KUTTA 

560 TFST=0 
555 DO 24 1=1 ,N 

24 SI VAR ( I ) = VAR ( I + 1 ) 

575 CUVAR ( 1 ) = V AR ( l ) 

576 DO 25 I = 1 ,N 
SDY ( I ) =DFR ( I +1 ) 

25 CUVAR <I + l)=SlVAR<I) + ( DX*0ER ( I +1 > )/2 • O 
CUVAR ( 1 ) =CUVAR ( 1 ) + DX/2 • 0 

CALL DFRSUR 

IF (II ,PQ. 4 ) GO TO 120 

580 DO 26 I = 1 * N 

SDY ( I ) =SDY (I ) +2 • 0*DFR ( I +1 ) 

26 CUVAR ( I + 1 ) =S 1 VAR ( I ) + ( DX*D5R ( I + 1 ) )/2 • 0 
CALL HFRSUB 

IF (II .FQ. 4 ) GO TO 120 

585 DO 27 I = 1 ,N 

SDY ( I ) = SDY ( I )+2» 0*DKR ( I +1 ) 

27 CUVAR ( 1 + 1 )=S1VAR< I )+DX*D£R( 1 + 1 ) 

CUVAR ( 1 )=CUVARU ) +-DX/2 • 0 

CALL D(=RGUB 

IF (II .^O. 4 ) GO TC 1 2G 

590 DO 90 I = 1 « N 

SDY ( I ) = ( SDY ( I )+n;R( I + 1 ) ) /6.0 
90 CONTINUF 

IF (LOOP) 28,28*29 

28 D030 I = 1 ,N 
SDY I ( I ) =SDY ( I ) 

Y I NCR ( I ) =0 .0 

30 DER ( I + 1 ) =FDERV ( I ) 

DX = DX/2 *0 

L OOP = 1 
GO TO 575 

LOOP WAS NOT 7FR0 

29 DO 31 I = 1 ♦ N 

31 Y I NCR ( I ) =Y I NCR ( I ) + SDY ( I ) 

IF (LOOP .EO* 2) GO TO 33 
DO 32 1=1 fN 

SI VAR ( I ) = V AR (1 + 1 ) +DX*Y I NCR ( l ) 

32 CUVAR ( I + 1 ) =S 1 VAR ( I ) 

CUVAR ( 1 )=VAR(1 ) 4-DX 
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LOOP=2 
CALL HPRSUB 

IF (IT .^O. 4 ) CO TO 120 
GO TO ^76 

33 LOOP=0 
H=?.0*OX 

FO 34 I = 1 iN 

FR VO VH ( I )= (YINCP( I )/2.n-S0Yl ( I ) )/lB.O 
ERRVAL ( I )=H#cRVnVH ( I ) 

34 51 VAR ( I ) = 5 1 VAR ( I )+DX*SDY ( I ) +FRRVAL ( I ) 

C 

C 51 VAR HOLO THF APPROXIMATE AN8WFRS 

C 

IF ( 5 C I M A X ) 

36 IF CAR5(5CIMAX-CI ) .LT, 1 *OE-12)GO TO 3R 
I F ( ARS ( H-C. I > .GT * 1 .n~- 1 p )GO TO 38 
or on^=o 
GO TO 60F 
38 DCODF=l 
605 CONTINUE 
I =0 

40 1=1+1 

IF (I .GT. N ) GO TO 45 

IF ( ABS ( S 1 VAR ( I ) ) .LT. EL£2 < I) ) GO TO 40 
R F L F R = A F 8 (FRRVAL ( I ) /^ 1 V AR ( I ) ) 

IF (RFLER mGTm cFLFKI)) GO TO 615 
IP (RFLFR .GT. RFL^IN(I)) OCOOF=1 
GO TO 40 
45 CONTINUF 

IF (DCOOF-1 ) 610,620,610 
610 CONTINUE 

IF ( S SPEC ) 4 1 ,42,41 
42 IF (SCI MAX >41 *43*41 
A3 CI=2.0*CI 
TEX (8 > = 1 
NHAF=NHAF-1 
GO TO 620 

41 IF (2* 0*ABS ( C I ) .LF. ABS(SCIMAX)) GO TO 43 
44 C I = 5 C I M A X 

T P X (8 ) = 1 
GO TO 620 
C 

C HALF INTERVAL 

6 1 R NHAf-" = MHAF+l 

TFX ( 9 ) = 1 
NV AR= I + 1 

J F (NHAF-8 ) 47, 47*505 
A 7 IF (LTFH .FQ, 0) GO TO 48 
TF?.T=I 

LT5H = o 

nflt=nflt- 1 

N- LTL=N-LTL+ 1 

FL TK = K 1 #S C LT (NFl T ) 

rfma I N=0 *0 

48 IF (T^SH • c O • 0) GO TO 49 
T F 5 T = 1 

TP = TP— SSPFC 
TPK = K 1 *T=> 

TPSH=0 

49 IF (SSPFC. O.n) GO TO 998 

TF ST = n 
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IF ( ABS ( C I -2 .n-H-DX ) • G T • 1 • F— 1 2 ) GO TO 1100 

998 C I =DX 

999 DX = DX/2 . O 
C IK=K1 *C I 
DO SO I = 1 ♦ N 

SI VAR( I ) =V AR ( 1 + 1 ) 

DER ( 1 + 1 )=FDFRV ( I ) 

SDY1 ( I )=YINCR( I ) — SOY < I ) 

50 Y I NCR ( I )=n.O 
KK3 = 2 

IF (I TEXT • FQ • 1) GO TO 637 

99 LOOP = 1 

GO TO 575 

1 1 00 CONTINUF 

IF (NHAF • GT • 1 ) GO JO 999 

NTS=NTS+1 

IF (NTS .GT. 13) GO TO 998 
ACV = VAR (1 )+C I 
ACVK=ACV*K 1 

IF (NELTL .EQ. O) GO TO 1102 
NLT = NELT 

1101 FLTK 1 =SELT (NLT ) *K 1 

IF ( AC VK .LT. FLTK1 ) GO TO llO^ 

NL T = NLT + 1 

IF (NLT .^Q. NTC) GO TO 1106 
GO TO 1 1 O l 
1 1 OP SFLT (NFLT ) =ACV 

00 TO 1 1 OS 

1 1 03 NL TP 1 = NL T + 1 

1 =NTS 

1 1 08 SFLT ( I ) = SFL T ( 1-1 ) 

I F ( I .FQ. NL TP 1 ) GO TO 1106 
1 = 1-1 

GO TO 1 1 08 
1 1 06 SF{_t (NLT ) = AC V 

1 1 OF nfltl = nfltl+ 1 
TFX (9 ) =0 
TFX ( 1 O ) = l 

FLTK=K1 * SFL T (N^l.T ) 

GO TO 999 
O 

C DOUBLE PRFCISIION UPDATING 

c 

620 LOOP = 0 
DH = H 

DO SI I = 1 ,N 

PHI=FRVOVH( I ) + Y I NCR ( T )/2.0 
DP H I =PH I 

51 CUVAR ( I + 1 ) = VAR ( T + 1 )+DH*DPH I 
CUV AR ( 1 )=VAR ( 1 ) +DH 

CALL DFRSUB 

IF (II .PQ. 4) GO TO 1 20 
CALL CHSUR 
IF ( I 1-2 ) 54 * 60 0 . 121 

121 TFST = 0 
S4 OO 57 1 = 1 « N 

S7 FDERV ( I ) = DER ( 1 + 1 ) 

SUMHAF = SUMHAF + NHAF-5TFP { 1 ) 
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ST FP ( 1 ) =STFP (P ) 

STEP ( P )=STEP ( 3 ) 

ST*~“P ( 3 ) =NHAF 
NHAF=0 
I FRR= 1 

IF (SUMHAr-8) 63,63,510 
6 3 DO F9 I = 1 , NIP 1 
69 VAR ( ! ) =CUVAR ( I ) 

TFX'( 1 2 ) = 1 
501 KK3=4 

IF (ITFXT *fq. \) GO TO 637 
68 IF (TFST .EO. } ) GO TO 530 
130 RETURN 


c 

C RFCOMPUT c interval 

c 

6 or tf ST = 0 

K'HA^rO 
11=1 
DX = C I 
TF X (1 1 )=1 
K<3 = 3 

IF (ITFXT • F Q # i ) GO TO 636 
7 n ClK = ri*ISl 

HO 60 I = 1 ,N 

D F R ( I + 1 ) =FDFRv ( T ) 


6~ CUVAR ( I ) = VAR ( I ) 

CUVAR(N+1 )= V A P ( N + 1 ) 

IF ( TPFH .EO. O) GO TO 61 
TP=TP— SPEC 
TPK=TP*K 1 
TPSH=0 

61 if (LTSH .EG. O) go to 655 
N^LT = NF[_t~ 1 
RF(V»A I N = 0 • O 
N**LTL=NFLTL + 1 

TK = SFLT fM c LT ) *K 1 
GO TO 5FK 


636 
635 

637 


1 OJ? 
330 


PRINT 1 83, VAR ( 1 ) , DX 
GO TO 103 

IF ( T F X ( 1 ) .EO. 1 )PRINT 
IF (TFX(3 ) .EG. 1 ) PR I NT 
I cr (TFX(3) *FQ. 1 ) PRINT 
IF(TEX(4).EQ.l ) PRINT 
I F ( TFX ( 5 ) •EQ. 1 > PRINT 
IF (TEX (6 ) .EQ. 1 ) PR I N T 
I c ( TFX(7).FG.l ) PR I NT 
I F (TFX ( 8 ) .EG. 1 ) PR I NT 
I F (TEX (9 ) .EG. 1 ) PR I NT 
I ^ (TFX ( 1 O ) .“G. 1 ) PR I NT 
IF (TEX( 1 1 ) .FO. I ) PR I NT 
IF (TFX( 12) .FO. 1 ) PR I NT 
I F ( TE X ( 13) •FO. 1 ) PR I NT 
IF ( T t X ( 14) .EG. 1 ) PR I NT 
I F ( TFX ( 1 5 ) .EG. 1 ) PR I NT 
DO 330 1=3.1 3 

TFX ( I ) =0 


171 4 V A R ( 1 ) 

172«CI *0 1 MAX* SPEC 
1 73 
1 74 . H 
1 754H 

1 76 4 H 
184, H 
177, Cl 
1 78 1 N V AR , C I 
1 86 , N V AR , OX 
1 83, VAR (1 ) , DX 
1 79 , VAR ( 1 ) 

1 80 
1 81 
1 8? 
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GO TO ( IPO ♦ng, 70, *=58 ) ,KKT3 

171 FORMAT C33H INITIALIZATION STARTS AT V A R ( 1 ) = « F 1 6 • R / ) 

172 FORMAT (4H C I = , F 1 5 • 6 , OH C I M A X = » r lF,RiBH S P F C = »F 1 F *8/ > 

173 FORMAT C37H DX TS T HF FULL COMPUTING INTERVAL C.I/) 


1 74 

format 

1AL VALi 

(28H 
JE/ ) 

DX 

I S A 

SHORTENED 

I NTFRVAL 

♦E15.8*25H 

DUF 

TO 

A 

CRITIC 

1 75 

FORMAT 
1 A LUF / ) 

(28H 

DX 

IS A 

SHORTENED 

I NTERVAL 

f E15*8<21H 

DUE 

TO 

A 

SPEC 

V 

1 76 

FORMAT 

(28H 

DX 

IS A 

SHORTENED 

I NTFRVAL 

,Ei5.Bt39H 

DUF 

TO 

BOTH A 

s 


lP^C AND CRITICAL VALUF/ ) 

177 FORMAT (27H Cl HAS Bffn LENGTHENED TO iF16*8/) 

173 FORMAT (5H VAR(,I2,32H) HAS CAUSFD Cl TO RE HALVFO TO ,E16.8/) 

179 FORMAT (2 ?H VAR ( 1 ) HAS BFEN UPDATFO TO ♦ F 1 6 • 8 « / ) 

180 FORMAT ( 3 1 H FRROR RFTURN-ELT NOT MONOTONIC/) 

181 FORMAT (55H FRROR RETURN-HAVF HALVED 9 TIMES OVER LAST 3 INTERVALS 
1 / ) 

182 FORMAT (45H frroR R-TURN-HAVF HALVED 9 CONSECUTIVE TIMES/) 

183 FORMAT ( 3 1 H I NTFRVAL RECOMPUTED AT VAR ( 1 )=iE16*8»9H WITH DX=«F16#B 
1 / ) 

184 FORMAT ( 25H DX IS SHORTENED I NTERVAL < E 1 6 • 8 , 28H DUF TO A PREVIOUS FL 
IT VALUF/) 

185 FORMAT ( 5H VAR(,I2«32H) HAS CAUSED DX TO BE HALVED TO «E16.8,38H Bu 
IT NOT Cl S I NCF Cl ALREADY SHORTFNF D/ ) 

500 IFRR=? 

TEX ( 1 3 ) = 1 
TF ST = 0 
GO TO 63 

505 I frr = 3 

TFX ( 1 5 ) = 1 
TF ST = O 
GO TO 501 

510 I FRR = 4 
TF ST = 0 
TFX ( 1 4 ) = 1 
GO TO 63 
END 
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Sample Output 

The following is an example case with the input quantities and results printed 

SNA Ml 


A 

= 

o 

• 

o 

SIGMA 

— 

0.1E-01, 

ALPHA 

= 

C* 2F+CC, 

PRMIN 

= 

0.5E-QU 

CII 

= 

0.2E-12, 

C IMAX 

= 

C. 4E-C3 , 

ALL U I 

= 

0.1E-04, 

SUL U I 

= 

C.12E+C1, 

TELTX 

= 

0.1E+00, 

EPI 

= 

0*1 E-0 6 , 

EP2 

= 

0. 1F-C6, 

NT 

= 

1* 

I CP BET A 

= 

0, 

$END 




$NAM2 




VARO 

il 

o 

• 

o 

* 

o 

• 

C, 0. SS5E- 

05, -0.199E-C4, 

ELE 1 

= C. IE-06, 

0.16-06, 

0.1 E-C6, 

ELE2 

* 

r- 

0 

1 

UJ 

T—i 

• 

o 

II 

0.1E-C7, 

0.2E-C7, 

ELT 

= C.1E+C1, 

I , I , I , I 

« It It It It I, 

SEND 





SNA M3 

SETA = 1, 
SENC 
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[ 


BODIES OF RE VC LI 7 I C N HAVING MINIMUM VISCCUS PLUS 


(NEWTONIAN) PRESSURE DRAG IN HYPERSONIC FLOW 


INPUT 

VO 


A 

SIGMA 

ALPHA 


9.95000CCCE-C6 

0 . 


1.00000000E-02 

2. 00000000 E-01 


INITIAL U OBTAINED FROM EULER EQUATION 
BETA NOT INCLLOED IN LAMBDA DOT EQUATION 


X 


Y 


Y PRIME 


0 . 

5* 0046 53 75E-C2 
1 *00 44653 8E- 01 
1* 5C 64653 8E-01 
2*01 246 5 38E-C 1 
2* 51 6465 38E-01 
3*02046 5 38E-01 
3. 52446538E-C1 
A. 028465 38E-C1 
4* 53246 538E-01 
5* C 3 6465 2 8E- C 1 
5* 54046538 E-01 
6* C4446 5 38E- C 1 
6*548465 38 £-0 1 
7* C5 246 5 38E-0 1 
7* 55646538 E -01 
8.C6C46538E-C1 
8*564465 38 E-0 1 
5* 066465 3 BE- C 1 
9*57 2465 38 E-0 1 
1*COCCCOOOE+O0 


9.950C0000E-C6 
5 • 2 4 6 885 54 E-0 3 
8.83052659E-C2 
1 • 19701598 E-02 

1 * 4 6 5 1 54 5 IE- C 2 
1.75564510 E-02 

2 * 0 1 2 7 6 5 4 7E- C 2 
? *2 59 309 39 E-0 2 
2.45714527E-C2 
2*7 2763376 E-0 2 
2. 951776 3 2E- C 2 

3 *17036372 E-0 2 
3.364C2272E-02 
3.59326325E-C2 
3. 7565C745E-C2 

4 *0001 1000 E-02 
4.198372 62 E-02 
4.35355483E-C2 
4* 5 6 5 8 81 5 7 E-02 
4*77 55513 IE- C 2 
4 • 93449 1 1 6E-C2 


5* 5599998 6E - 01 
7.82930417E-02 
6. 574G673 IE -C 2 
5 *93708 350 E-0 2 
5. 52342547E-C2 
5.22272838E-02 
4* 5 89 34 04 4E -02 
4 • 800 231 59E-0 2 
4. 6422602 IE -02 
4. 50 72670 IE -0 2 
4.38985957E-02 
4*28630186E-02 
4. 15390632 E-02 
4. 11067857E-C2 
4. C351 01 07E-02 
3*9659942 5E -02 
3* 90242 53 CE-02 
3* 643643 80E- 02 
3 .7890364 1 E-02 
3. 7 3809887E-C2 
3 *69745 149E-02 


TOTAL DRAG COEFF 
4*C5 765800 E -03 


PRES DRAG COEFF FRICTION DRAG COEFF 
4 *09 76 5 800 E-0 3 0. 


GO 
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Figure 1 - Drag coefficient ratio of minimum-drag bodies with given length and volume, o = 0.01. 
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Figure 8.- Minimum-drag body profiles. Length and volume are given 
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